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PACS 45 . 70 . -n - Classical mechanics of granular systems 
PACS 82 . 70 . Rr - Foams 

Abstract. - We experimentally probe the vicinity of the jamming point J, located at a density 
(j) corresponding to random close packing {(j)rcp — 0.842), in two dimensional, bidisperse packings 
of foam bubbles. We vary the density of the foam layer and extract geometrical measures by 
image analysis. We confirm the predicted scaling of the average contact number Z with and 
compare the distribution of local contact numbers to a simple model. We further establish that 
the distribution of areas p{A) strongly depends on cj). Finally, we find that the distribution of 
contact forces p{f) systematically varies with density. 



Since the seminal work by Bolton and Weaire [T soft 
frictionless discs or spheres have become the Ising model 
for the Jamming transition pl-[4]. Jamming is believed 
to capture the transition between flow and arrest in a 
wide variety of disordered media such as foams, emulsions, 
granular media and (colloidal) suspensions, and this idea 
has lead to an upsurge of simulations revealing critical be- 
haviour as a function of the distance to the critical point 
J, which is reached precisely when the applied pressure 
vanishes [sIIsIIg]. While some of these predictions have 



been tested experimentally |7- 10 , many others still await 
experimental verification. 

While the original incarnations of the soft sphere model 
explicitly make the link to foams [ip] , this connection has 
not been explored in recent times. Nevertheless, foam bub- 
bles, as well as emulsion droplets are the closest physical 
analogue of frictionless spheres. Their elastic interaction 
is close to that of a linear spring pT|p!2] and solid friction 
is absent. Under flow, only a velocity dependent viscous 
friction, that is now fairly well understood [5]|6[[l3J[l4] acts 
on the bubbles. 

In this Letter we experimentally probe the behaviour of 
jammed packings near point J. We do this by generating 
two-dimensional packings of foam bubbles. We vary the 
packing fraction, (/), and for each density generate many 
distinct static packings of foam bubbles. We then extract 
various statistical and topological quantities through im- 
age analysis and investigate whether these measures signal 



an approach of the jamming transition at (j)c as the density 
is varied. 

First, we investigate the scaling of the average contact 
number per bubble Z with the distance to (pc Our cen- 
tral result is that the contact number scales like 7j — T^c ^ 
{(j) — 0c)^'^- To this end, we first resolve the apparent dis- 
crepancy between our findings and simulations by point- 
ing out the differing ways of measuring cj) between simu- 
lation and experiment. We find (j)c to be located around 
(j)c = 0.84, in excellent agreement with previous predic- 
tions [l][3J|5J[T5] , and obtain, for the first time, a quantita- 
tive experimental confirmation of critical scaling at point 
J in an appropriate experimental system. 

Besides this global measure, we also investigate the dis- 
tribution of contact numbers at the bubble scale as a func- 
tion of global Z. We compare to a recent model |16 and 
find excellent agreement with no adjustable parameters, 
which is remarkable given that this theory was developed 
for frictional packings. 

The distribution of free available area per bubble, p{A) 
plays an important role in statistical mechanics descrip- 
tions of jammed matter [l7 19 . We show that the correct 
way of extracting p{A) is by tesselating the foam packing 
with the navigation map [2Q{[2T] . We fit the obtained dis- 
tributions of free area per bubble to gamma distributions 
18,22 , and show that, in contrast to granular packs, this 
distribution is far from universal: as we decrease the pack- 
ing fraction (j) towards point J, we find an excess of large 
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Fig. 1: (Color online) Top (a) and side (b) view of the setup. 
(a) Bubbles are contained in a circular reservoir and covered by 
a glass plate. Images are taken in the highlighted area, (b) A 
stepper motor drives a rod that stirs the bulk solution, thereby 
rearranging the packing, (c) From left to right: raw image, raw 
image with reconstructed bubble areas, reconstructed bubble 
areas, from which (j) is calculated, and contact network, from 
which Z is found. Scalebar denotes 5 mm. 



available areas, and in conjunction with this an increase 
of the compactivity. 

Finally, our navigation map tiling allows us to extract 
the deformed interfaces between adjacent bubbles. These 
facets are proportional to the contact forces, and we thus 
extract the distribution of the contact forces p(/). We find 
that its tail is close to exponential for (j) close to (/)c, and 
becomes steeper for large 
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Setup — We prepare a surfactant solution consisting of 
0.5 % volume fraction Dawn dishwashing liquid and 15 
% glycerol in demineralized water (viscosity rj = 1.8 ± 0.1 
mPa-s and surface tension cr = 28±1 mN/m) in a large cir- 
cular reservoir (r = 190 mm) of depth 30 mm, see Fig. Ilia). 
A bidisperse (50:50 number ratio) bubble monolayer is 
produced by flowing nitrogen through two syringe nee- 
dles immersed at fixed depth in the soapy solution. The 
resulting bubbles of 1.8 ± 0.1 and 2.7 ± 0.1 mm diam- 
eter are gently mixed to produce a disordered bidisperse 
monolayer and are covered with a 10 mm thick glass plate, 
see Fig. Ilia). The weighted average bubble diameter (d) 
is 2.25 mm. We light the bubbles slanted from below and 
cover the bottom of the reservoir with a black plate, to 
enhance contrast. 

The bubbles contact the top plate, see Fig. [ijb), which 
is completely wetted by the soap solution, and the liquid 
fraction of the foam can be varied by varying the distance 
between glass plate and liquid surface between 3 and 0.2 
mm. This in itself is not a proper measure of (j) since the 
relation between (j) and the gap is strongly hysteretic — 
(j) depends not only on the gap, but also on an uncon- 
trolled confining pressure. Therefore, we will determine (j) 
from experimental images, in a procedure outlined below, 
see Fig. [Tic). We check that coalescence, segregation and 
coarsening are negligible. 

Experimental protocol — A stepper motor is glued to the 
glass plate and is connected to an aluminum rod through 
a hole in the glass plate, see Fig. [ijb). We agitate the 



surfactant solution underneath the bubbles by driving the 
rod back and forth with an amplitude of 1 radian, its angu- 
lar velocity alternating between +0.6 and -0.6 radians per 
second. We emphasize that while agitating the bulk solu- 
tion leads to strong mixing of the packing, we inject little 
enough energy to avoid bubble break-off, as evidenced by 
the absence of satellites. After 4 oscillations we stop the 
motor, after which the packings slowly relax to a mechan- 
ically stable state. We probe the relaxation of packings at 
varying (j) to determine the waiting time between agitation 
and image acquisition. To this end, we record sequences 
of images with a CCD camera and measure the variance of 
the intensity fluctuations of all pixels in difference images. 

After waiting for this time, which is of the order of min- 
utes, we record one image in the previously agitated region 
with a 6 megapixel photo-camera (Canon 20D). The im- 
age contains between 350 and 700 bubbles, depending on 
the packing fraction cj). For various fixed gaps between liq- 
uid surface and glass plate we repeat this procedure 100 
times and thus obtain 100 packings at roughly equal pack- 
ing fraction. We visually inspect the resulting packings to 
be distinct in appearance, and it is these images that we 
analyse in the following. 

Determining (j) — We extract our crucial control param- 
eter (j) from the experimental images by advanced image 
analysis, see Fig. [2|c): we first binarize the image, after 
which both the bubble centers and the interstices appear 
bright. We then remove the interstices by morpholog ical 
operations and end up with an image consisting of bright 
bubble centers. We dilate these centers and add up a neg- 
ative of the original binary image — in which the bubbles 
appear as bright rings — to arrive at the final image, in 
which the bubbles are represented by bight discs against 
a black background. From this image we can readily cal- 
culate the area fraction (j). Note that, in principle, the 
concept of packing fraction is problematic for a monolayer 
of three-dimensional bubbles. We choose our lighting of 
the bubbles such that the contacts between adjacent bub- 
bles are optimally resolved. In other words, we image a 
slice from the packing where the bubbles are the broadest 
and calculate a 2D packing fraction from this slice. 

Scaling of Z with A(j) — We first determine the scaling 
of the average contact number Z with 0. To determine 
Z, we locate the center of mass of each bubble in the im- 
age, and after Delaunay triangulation and a subsequent 
removal of bond vectors for non-touching bubbles, we ob- 
tain the contact network of the bubbles in the image, from 
which we calculate Z, see Fig. [Tic). 

Our results are presented Fig. [2J where we plot the val- 
ues for each distinct packing (grey dots), the average over 
all 100 images for each packing fraction (black circles) as 
well as a powerlaw fit of the form Z = 4 + Zo(0 — ^c)^ (red 
solid line), where 4 is the contact number at isostaticity. 
The best fit gives us Zq = 4.02 ± 0.20, (/)c = 0.842 ± 0.002 
and P = 0.50 ± 0.02, in remarkable agreement with the- 
oretical predictions by O'Hern et al and Durian. [2j|3] 
who found Zq = 3.6 ± 0.5, ^c = 0.841 ± 0.002 and 
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Fig. 2: (Color online) Contact number Z of packings versus 
their packing fraction cj). Grey scatter: data for every indi- 
vidual image. Circles: data averaged over experimental run 
at approximately constant packing fraction. Solid line is fit to 
Z = 4 + Zo((/)-0c)^, withZo = 4.02 ±0.20, 0c = 0.842 ±0.002 
and /3 = 0.50 ± 0.02. Upper inset: same data on log- log scale. 
Lower inset: Z versus experimentally determined packing frac- 
tion (j)exp- The fit has a power law exponent of 0.70. 



/3 = 0.49 ± 0.03. 

Note however, that the range of packing fractions we can 
scan over, extends to a surprisingly large value of = 1.06. 
This is due to a striking discrepancy between the manner 
in which (j) is calculated in simulations and in experiments. 
In simulations, the area or volume of spheres is fixed, and if 
one knows the number of particles in the periodic box, (j) is 
readily calculated. In experiments, however, (j) can only be 
inferred from experimental images. This difference results 
in the following: if particles overlap, the overlapping area 
of the two particles is counted twice in simulations, while 
it is only counted once in our experiment. This doubly 
counted area scales with the overlap ^ as Aqv ~ ^^^^, which 
stems from the fact that the deformed area scales as Vc x 
i = VT X i 11.- Since i ^ {(j) — ^c) 4 , the conversion 
between a packing fraction extracted from a simulation (pth 
and its experimentally accessible counterpart (pexp should 
read: 



(t>exp = (t>th - C{(j)th - </>c)^^^. 



(1) 



We calculate both 



and 



^exp <^^^^ 4^th from numerically gener- 
ated packings, and determine the pre- factor C = 0.95. We 
then invert Eq. ([l]) and calculate the (l)th corresponding to 

our (t)exp' 

When plotting our data against (j)th as in Fig. |2J we 
excellently match simulations, while we find an apparent 
scaling exponent j3 = 0.70 if we plot Z as a function of the 
experimentally determined (l)expi see lower inset of Fig. l2| 
owing to the nontrivial relation in A0 between 0^^ and 

y^exp- 



We are not the the first to experimentally investigate the 
scaling of Z with (j). Majmudar et al. 9 have extracted 
the same quantities from images of two-dimensional, fric- 
tional, photoelastic discs and compared these to predic- 
tions from simulations. From their data it appears the 
prefactor Zq ~ 16, inconsistent with simulations. Our 
results do allow for a direct comparison with frictionless 
jamming predictions, which can be seen from the excellent 
agreement between parameters. 

Local contact fractions — Besides the average contact 
number per packing Z we can also extract the fraction Xz 
of bubbles in each image that has z contacts. We aver- 
age these fractions over all images that correspond to a 
global packing fraction (and contact number Z) cf. the 
black circles in Fig. [2] We plot these fractions versus the 
average Z in Fig. [3J we see clear trends in the abundance 
of contacts at the particle level, to which we apply a very 
recent model 16 . 

This model predicts the fractions of 4 species 
{xn^ '•'^Xn-\-3} in a packing, given the global Z and the 
variance a^ = X^iLn ^^(^ ~ 0^- This constraint, together 



with the trivial normalization constraints ^^ 



n+3 



1, 



E- 



n+3 



tXi 



Z and the ill-understood, but empirically ob- 



served [^ constraint that the number of particles with odd 
and even contacts is equal, leads to a set of of equations, 
the solution of which is: 



^n 


= ((Z-(n + 2))2 + a2-l/2)/4 


(2) 


Xn-\-l - 


= (-(Z-(n+l))2-(T2 + 5/2)/4 


(3) 


Xn+2 = 


= (-(Z-(n + 2))2-cT2 + 5/2)/4 


(4) 


^n+3 = 


= ((Z-(n + l))2 + a2-l/2)/4 


(5) 



Since we know Z and cr^(= 0.75) from the data we 
can obtain the fractions Xi without any free parameters. 
However, we measure non-negligible fractions of not 4, 
but 5 species. We therefore apply the model for n = 3 
to 4 < Z < 4.75, where X7 ^ and for n = 4 to 
4.97 < Z < 6 where X3 ?^ 0. We obtain good agree- 
ment between data and theory, see Fig. |3) for both ranges 
of validity. At high Z, the bidisperse nature of our pack- 
ing is visible: we observe an excess of both particles with 
5 contacts and with 7, while particles with 6 contacs are 
underrepresented. We have observed that large bubbles 
carry the majority of 7's, while small bubbles mostly have 
5 contacts at these values of Z (data not shown). This is 
natural in bidisperse packings, as this occurs whenever a 
large and a small bubble contact. It thus indicates an ab- 
sence of crystalline order, which would lead to an increase 
of 6's. Area distributions — We now turn to tesselations of 
our foam packings. These tesselations yield two important 
"connectors" between local geometry and global response. 
Firstly one can readily extract the distribution of avail- 
able area per bubble p{A) 19,24 28 , which serves as the 
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Fig. 3: (Color online) Fractions of bubbles in the foam with n 
contacts as a function of Z. Solid lines: solutions to Eqs. (2-5) 
for the species listed at the top of the graph. 



materials. Secondly, the size of the contact area's between 
bubbles can be measured, from which the distribution of 
forces p{f) can be extracted which is needed as an input 
into any prediction of the mechanical properties of pack- 
ings. 

The thermodynamical description of granular materials, 
as introduced by Edwards and Oakeshott 29 translates 
the concepts underpinning equilibrium thermodynamics 
to conglomerates of a- thermal particles such as bubbles. 
The volume V (in 2D the area A) takes on the role of 
energy, while a compact ivity x replaces temperature. 

Aste and Di Matteo have derived the form of the distri- 
bution function p{A^ k) with such an approach [l8]: 



p{A,k) 



k" 



[A 



k-\ 



{k-l)\{{A)-A^,^Y 



exp 



{A) - A„ 



(6) 

with k a shape parameter — which has been found to 

take on one universal value in granular packings, namely 

k = 12 [24!, Arain the minimum available area per bubble 

and {A) the mean. 

For deformable particles, it is not immediately obvious 
how the areas A should be computed. Here we calculate 
p{A) from our experimental images for each packing frac- 
tion both from a Voronoi and a navigation map tesselation, 
which are detailed below. Upon fitting Eq. (|6| to our data 
we obtain excellent matches to all distributions measured 
in both ways, but we find strongly nonuniversal behaviour 
in the parameter k for the navigation map method which 
we argue to be the correct one. 

We measure the probability distribution of free areas 
^(^4) by calculating the Voronoi area distribution of the 
grid of points that represent the centers of mass of the 
bubbles. However, Voronoi cell edges do in general not 
respect the bubble perimeter, see Fig. |4ja) and thus the 



Voronoi cell does not represent the free area per bubble. 
For hard spherical objects one can get around this problem 
by weighting the grid points according to the sphere ra- 
dius (Voronoi-Laguerre tessellation) [19^, however, in our 
experiment, the bubbles are not only bidisperse, but in 
general also deformed and the flattened contacts can be 
curved. 

To fully take the effects of both deformations and bidis- 
persity into account, we calculate what is called the navi- 
gation map 20,21 . In this method, we assign the intersti- 
tial area to those bubbles whose perimeter is closest. The 
result is shown in Fig.[4Fb): we obtain free areas per bub- 
ble that respect the bubble edges and follow the curvature 
of the contacts. 

With both methods we obtain bimodel distributions for 
A, which we split according to the size of the bubbles. Dis- 
tributions for the larger bubbles are shown in Fig. [4Fc,d) 
for Voronoi and navigation map tesselations respectively. 
Distributions for the smaller bubbles behave the same. 
For the Voronoi tesselation we find that the shape of the 
distributions is roughly independent of the packing frac- 
tion, with all distributions being slightly skewed (see inset 
of Fig. mc)), while for the navigation map tesselation we 
find that an excess of large available area develops near 
jamming, leading to strongly asymmetric distributions. 

We quantify these observations by fitting Eq. (|6| to our 
experimental distributions. We treat /c as a fit parameter 
and extract Amin from the data: it is the minimal hexagon 
that can enclose a bubble at a given packing fraction. Re- 
sults are plotted in Fig. |4je,f): the nearly constant shape 
of the Voronoi distribution is reflected by the near con- 
stant value of k we extract, with < k >= 13.1 for large 
bubbles and < k >= 11.4 for small bubbles when aver- 
aged over all values of (j). On the other hand, we observe 
a systematic trend in k for the navigation map distribu- 
tions, the nearly Gaussian distributions at high (j) can be 
fit when /c ~ 50, but k systematically decreases with de- 
creasing packing fraction (/) to a value of 6 near jamming. 

We thus observe that p{A) strongly depends on the par- 
ticular tesselation we use to extract it. For Voronoi tesse- 
lations we find a k which is remarkably close to the value 
found by Aste and coworkers in hard-sphere packings, 
while for the navigation map we find a strong variation 
of k with 0, that reflects the increasing amount of excess 
available area per bubble with decreasing packing frac- 
tion, which is also reflected by an increasing compactivity 
X = {{A) - Amin)/k iW, see inset of Fig. ^f). Since 
the navigation map tesselation respects bubble edges, we 
believe this tesselation to be more physically appropriate, 
and we see the decrease of k and the associated broadening 
of the tail of p{A) as a signature of the approach of point 
J, or, equivalently, the recovery of hard-sphere behaviour. 

The force distribution p(f) — By construction, the navi- 
gation map bisects touching bubbles at their contact area. 
As a result, we can extract the radius re of the deformed 
contacts between bubbles. Since the contact force /12 be- 
tween two bubbles is linear in the deformation, it is related 
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Fig. 4: (Color online) (a) Experimental image with Voronoi tesselation of the centers of mass of the bubbles. Note the intersection 
of bubbles by the Voronoi cell edges, (b) Navigation map tesselation of same image: the cell edges do not intersect the bubbles, 
(c) Area distributions from Voronoi tesselation for 4 different packing fractions, as stated in legend. Solid lines are fits to 
Gamma distribution Eq. p] Little trend with is seen, (d) Area distributions from Navigation map tesselation for the same 
packing fractions. Solid lines are fits to Gamma distribution Eq. p] Near (/)j, broad tails develop, (e) k extracted from fits to 
Voronoi cell distributions. No trend with can be seen, (f) k extracted from fits to Navigation map distributions. A strong 
variation of k with is visible. The inset shows the compact ivity x — ((^) ~ Amin)lk^ which increases towards point J as /c 
decreases. 



to Tn as 



30 



A 



R1R2 



(7) 



with Ri the radius of bubble i. We can thus extract the 
force distribution p{f) of interbubble contact forces. An 
experimental image with the contacting facets in white 
and the extracted contact forces in blue is shown in 
Fig. da). 

The forces are oriented along the bond vector between 
bubbles and thus we can check if the forces balance on 
each bubble. To check this, we plot the relative error 
IISf|l/X]l|f|l i^ Fig.jsja). We see that the error is about 
10 % for the densest packings, and increases for packings 
that are closer to jamming. 

In Fig.lsjb) we plot, for various (/), the p{f) of all forces 
from all images at that packing fraction. We convert Tc to 
forces using the blob radii Ri we obtain after binarizing 
and removing the interstices. These radii are proportional 
to, but smaller than the actual bubble radii. We there- 
fore plot the force distributions in arbitrary units. The 
distributions are peaked and exhibit broad tails, reflecting 
the heterogeneity in the forces. Their shapes are similar to 
those found in grains and emulsions pQ[|3T , In the inset of 



Fig. [5|b) that the peak of the distribution shifts towards 
lower forces as we approach point J, in accordance with 
[321. 



To obtain more information on the shape of the tail of 
p(/), we now rescale each p{f) from an individual image by 
its mean and average these distributions for each (j). Just 
adding up the force distributions for each frame, as was 
done for Fig. [5|b), washes out any signature of Gaussian 
tails 32 . The result is shown in Fig.jsjc): we now observe 
changes in the shape of the tail away from jamming — 
with increasing compression the decay of the tail seems 
to become steeper than exponential. This is illustrated 
by the inset of Fig. Isjc), where we plot ln{p{f))/f which 
tends to a constant for exponential tails and decreases for 
faster than exponential decay. We thus provide further 
evidence that one can find steeper than exponential decay 
in p(/), although we do not have sufficient statistics to 
conclude that the tails cross over to Gaussian tails away 
from jamming. 

Conclusion. We have experimentally investigated the 
behavior of soft frictionless discs near the jamming point 
at zero stress and temperature, by sampling distinct pack- 
ings of foam bubbles. Experimental data of systems near 
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Fig. 5: (Color online) For the same as in Fig. p] we calculate (a) Probability distribution of relative error in force balance 
on each bubble: for packings closer to point J, the error increases. Inset: experimental image with force network overplotted: 
flattened contacts in white, and magnitude of the forces indicated by thickness of vector and darkness, (b) Force distributions, 
rescaled by the mean / of all configurations: the peak shifts to smaller forces as one approaches point J (see inset), (c) Force 
distributions, for each configuration we rescale its p{f) by the mean. The average now exhibits steeper than exponential decay. 



jamming are still scarce and we hope that the various mea- 
sures that we have extracted can help to decide between 
various competing theoretical descriptions of jammed mat- 
ter. 
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